function [lastmaxdi] =SuffNSolBellf(PI,M,vbarv,bet,gam) 
% Check if sufficient condition for discounting is met

N = max(size(vbarv));

Ms = (PI.*M)*gam;
XI = eye(N);

K=50000;   % can make larger if want

maxdisc = zeros(K,1);

for k=1:K

    MXI =  Ms*XI ;
    PXI =  bet*gam*PI*XI;
    Ixi = MXI*ones(N,1)> PXI*ones(N,1);
    
    XI = PXI;
    XI(Ixi,:) = MXI(Ixi,:);

    maxdisc(k)=max(XI*ones(N,1));
end
 
lastmaxdi = maxdisc(K);
if lastmaxdi<1
    'Satisfies discounting condition'
elseif lastmaxdi>1
    'Discounting NOT Satisfied -- change parameters'
end
    
figure(10)
plot(maxdisc)
title('Largest discount factor function of K')

